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INTRODUCTION 


This  final  progress  report  details  key  theoretical  and 
applied  advancements  in  the  field  of  Moment  Closure 
discovered  through  the  course  of  this  project. 

Specifically,  the  objective  of  this  work  was  to  study  the 
application  of  moment  closure  methods  under  cumulant 
truncation  to  the  transient  analysis  of  the  sortie 
generation  process.  The  work  that  was  originally  proposed 
included  1)  investigating  the  accuracy  of  moment  closure 
procedures,  2)  developing  procedural  improvements  to 
increase  accuracy,  and  3)  applying  these  to  the  sortie 
generation  process.  As  the  work  progressed,  it  became 
apparent  that  these  original  objectives  needed  to  be 
augmented  to  accomplish  the  proposed  work.  As  such, 
additional  research  objectives  of  this  project  include  4) 
creating  a  computationally  efficient  moment  closure  program 
and  5)  investigating  the  stability  of  moment  closure 
methods . 

As  a  result  of  this  work,  the  following  main 
accomplishments  have  been  achieved  1)  the  creation  of  a 
self-contained  Mathematica  package  with  users  manual  that 
provides  moment  closure  approximations  for  large  stochastic 
networks,  2)  the  development  of  analytical  procedures  to 
check  for  network  stability,  3)  the  development  of  an 
optimal  truncation  policy  based  on  the  maximal  order  of 
polynomial  intensity  (rate)  functions  for  the  network,  4)  a 
loose  correlation  measure  of  the  error  between  the  accuracy 
of  moment  closure  approximations  and  the  traffic  intensity 
of  the  network,  and  5)  the  application  of  moment  closure 
methods  to  the  large-scale  sortie  generation  process.  Each 
of  these  accomplishments  will  be  discussed  individually  in 
the  body  of  the  report .  There  have  been  three  research 
publications  developed  as  part  of  this  work,  two  invited 
presentation,  and  two  contributed  presentations.  At  the 
present  time,  a  paper  titled  ""  is  in  the  process  of  being 
re-submitted  to  Naval  Research  Logistics. 

RESEARCH  FINDINGS 
Creation  of  Moment  Closure  Program 

At  the  beginning  of  this  project,  the  Mathematica®  program 
initially  authored  by  the  PI  for  the  analysis  of  stochastic 
systems  was  deemed  to  be  inefficient  for  analyzing  large 
stochastic  systems,  such  as  the  sortie  generation  process. 


As  an  example,  a  10-node  version  of  the  sortie  generation 
model  described  in  Dietz [1]  would  take  several  hours  to 
enter  and  would  not  fully  evaluate  due  to  memory 
constraints.  This  same  model,  however,  can  be  entered  and 
evaluated  in  a  matter  of  minutes  using  the  streamlined  code 
that  was  developed  through  this  project. 

The  process  of  streamlining  the  moment  closure  code  was 
ongoing  throughout  this  project.  Significant  milestones 
include  the  completion  of  a  preliminary  streamlined  version 
in  April  2004,  the  completion  of  a  user-friendly  version  of 
this  code  in  August  2004,  and  the  ultimate  packaging  of 
this  code  in  March  2005.  This  program  provides  for  the 
efficient  analysis  of  arbitrarily  large  networks  and  only 
requires  user  to  input  the  intensity  functions  of  the 
network,  the  desired  truncation  level,  and  the  time  horizon 
over  which  the  model  is  to  be  evaluated.  The  contents  of 
the  Mathematica  package  may  be  found  in  Appendix  A  of  this 
report,  a  users  manual  in  Appendix  B,  and  a  sample 
implementation  in  Appendix  C.  This  Mathematica  code  may  be 
retyped  and  saved  as  a  .m  package  file,  or  an  electronic 
version  may  be  obtained  from  http://web.nmsu.edu/~-tmatis. 
The  majority  of  the  coding  work  was  performed  by  Karl  Adams 
and  Ivan  Guardiola  under  the  direction  of  the  PI,  and  the 
users  manual  was  co-authored  by  Amara  Nance  and  the  PI. 

Stability  of  Moment  Closure  Methods 

"Moment  Closure"  is  a  method  used  in  order  to  close  the 
open  set  of  Differential  Equations  obtained  from  the 
transition  intensities  using  the  Random  Variable  technique 
that  define  the  moments  of  the  system.  Currently  two  forms 
of  closure  are  being  implemented  throughout  the  scientific 
community.  These  two  forms  of  closure  are  neglect  and 
parametric.  Under  the  closure  method  of  "neglect"  we  draw 
to  the  assumption  that  all  moments  and  cumulants  above  a 
certain  "closure  level"  are  insignificant  in  the  sense  that 
they  are  not  vital  in  order  to  approximate  the  low  order 
cumulants  of  interest.  The  closure  method  of  "parametric" 
is  done  by  making  an  assumption  of  the  underlying 
statistical  distribution  and  formulating  the  high  order 
moments  and  cumulants  in  terms  of  the  lower  ones.  The 
underlying  statistical  distribution  gives  us  an  insight 
into  how  the  high  order  moments  should  approximated.  For 
example,  by  making  the  assumption  that  the  underlying 
distribution  is  "Poisson"  we  would  set  the  high  order 
moments  and  cumulants  above  n>2  to  be  equal  to  the  first 


cumulant,  this  is  due  to  the  characteristics  of  the  Poisson 
distribution.  These  parametric  assumptions  allow  us  to  use 
stability  analysis  to  fully  interpret  the  system's 
solutions  both  in  transient  as  well  as  in  steady  state. 

Stability  analysis  consists  of  the  deriving  the  system's 
stationary  points  otherwise  known  as  the  critical  points  or 
equilibrium  points  of  the  system.  This  is  done  by  solving 
the  system  of  differential  equations  after  the  parametric 
assumption  is  made  in  the  same  manner  as  solving  a  linear 
system  of  equations.  A  brief  mathematical  procedure  summary 
consists  of  finding  the  critical  points  of  the  system  of 
equations,  deriving  the  Jacobian  Matrix  at  each  of  the 
stationary  points,  deriving  the  Eigensystem,  which  contains 
both  the  Eigenvalues  and  Eigenvectors  of  the  system  at  each 
point,  and  deriving  a  staring  point  for  manifolds  that  will 
divide  the  phase  portrait  into  feasible  and  unfeasible 
regions.  The  procedure  above  can  be  done  by  choosing  the 
appropriate  closure  of  the  high  order  moments  and  cumulants 
based  on  an  underlying  statistical  distribution  assumption 
and  setting  the  differentials  to  zero.  This  will  yield  the 
corresponding  critical  points  of  the  system.  Stability 
analysis  continues  by  then  giving  a  description  and 
classification  of  the  system's  critical  points.  The 
Jacobian  of  the  system  will  then  give  us  a  means  of 
obtaining  such  classification  of  the  critical  points  as 
"source,"  "sink,"  or  "saddle".  Consider  the  following 
system 

t  =  f(x,y) 

t=g(*,y) 

Looking  at  this  system  asymptotically  we  set  the 
differentials  to  0.  Thus, 

0=f(x,y),0  =  g(x,y) 

We  then  find  the  equilibrium  or  critical  points  of  the 
system.  Suppose  that  (r0,y(,)is  an  equilibrium  point.  Then 
let, 

J=  i(xo’y0 )  f (wJ 
jsKwJ  %(xo>y0)_ 

be  the  Jacobian  Matrix  evaluated  at  the  equilibrium 
point  (x0,y0)  .  The  Eigenvalues  are  obtained  by  solving 
Det(J-AI)  =  0  in  terms  of  A.  The  Eigenvalues  allow  us  to 
fully  explore  and  classify  the  equilibrium  points.  These 
classifications  can  be  determined  if  the  Eigenvalues  follow 


the  following  criterion.  If  the  Eigenvalues  of  the  Jacobian 
matrix  are  negative  real  numbers  or  complex  numbers  with 
negative  real  parts,  then  the  equilibrium  point  of  the 
system  is  classified  as  a  stable  "sink"  or  "spiral  sink." 
This  solution  will  approach  this  point  as  t-)».  If  the 
Eigenvalues  are  positive  or  complex  with  positive  real 
parts  then  this  solution  will  move  away  from  this  point  as 
t— ><».  Thus,  such  a  point  is  an  unstable  point.  The  point 
can  then  be  classified  as  a  "source"  or  "spiral  source."  If 
the  Eigenvalues  are  positive  and  negative  parts  then  this 
point  is  classified  as  a  "saddle." 

Manifolds  are  solutions  of  the  system  of  differential 
equations  which  determine  the  behavior  of  solutions,  within 
all  feasible  space.  These  manifolds  are  of  importance  do  to 
the  fact  that  the  behavior  of  solution  on  either  side 
differ.  For  example,  on  one  side  of  the  manifold  solutions 
will  tend  to  a  stable  equilibrium  point,  where  as,  on  the 
other  side  solutions  might  tend  to  some  infinite  value. 
Therefore,  manifolds  will  be  employed  in  order  to  separate 
the  phase  portrait  into  feasible  and  unfeasible  regions. 

The  following  example  considers  a  single  node (compartment) 
Markov  system.  Let  the  immigration/birth  rate  function  be 
given  by  fx  (X)  =  A  =  4,  and  let  the  death  rate  function  be 
given  by  f-i  (X)  =  px2  =2x2.  The  system  is  graphically 
depicted  as 


We  find  the  following  results  for  the  stability  of  the 
above  system  under  a  normal  and  Poisson  parametric 
assumption,  yielding  the  following  closure  procedures. 

Normal  Distribution:  ki(t)=0  V  i>3 
Poisson  Distribution:  ki(t)=kj(t)  V  i,j>l 

Phase  plots  for  these  systems  are  given  in  Figures  1  and  2 . 

The  results  of  stability  analysis  in  moment  closure 
approximations  are  important  for  various  reasons.  First, 
through  this  analysis  we  are  able  to  determine  a  domain  for 
initial  conditions  in  which  solutions  will  behave  in 
accordance  to  the  birth  death  model  characteristics  and 
expectations.  This  domain  determines  which  initial 
conditions  will  yield  relevant  results  and  even  determine 
whether  solutions  can  be  attained  or  not.  Through  a  study 


Normal  Assumption 
x-axis=k1(t),  y-axis=k2(t) 


Equilibrium  Points: 

1.  {1.61803.-0.618034} 

2.  {1.1} 

3.  {0.618034,1.61803} 

Classification: 

1.  Unstable 

2.  Stable 

3.  Saddle 
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Figure  1 :  Stability  Analysis  under  a  Normal  Assumption 


Poisson  Assumption 
x-axis=k1(t),  y-axis=k2(t) 


Equilibrium  Points: 

1.  {-1.68614,-0.84307} 

2.  {1.18614,0.59307} 

3.  {0,2} 

Classification: 

1 .  Unstable  point.  “Source” 

2.  Stable.  “Sink” 

3.  Saddle 


Figure  2:  Stability  Analysis  under  a  Poisson  Assumption 


of  manifolds  we  are  able  to  graphically  interpret  such  a 
domain  as  shown  in  the  previous  figures.  The  manifolds  act 
as  separators  of  solutions  behavior  from  one  side  to 
another  we  can  clearly  see  that  the  solution  tend  toward 
very  different  values.  This  analysis  shows  us  another  mode 
of  establishing  error  involved  with  the  approximation 
method  based  on  making  certain  parametric  assumptions. 

This  is  an  important  aspect  that  has  not  been  fully 
explored  by  any  others  in  the  scientific  community.  Moment 
Closure  methods  are  currently  being  employed  in  turbulence 
air  flow  models,  signal  processing  and  many  other  in  which 
people  are  employing  the  moment  closure  method  specific  to 


parametric  assumptions.  These  parametric  assumptions  should 
be  under  full  exploration  and  determination  of  system 
behavior  prior  to  implementation  of  closure  under  such 
assumptions.  We  explore  this  issue  in  order  to  assure  that 
certain  assumption  carry  certain  error,  stability  and 
solution  behavior. 

Optimal  Truncation  Policies 

The  accuracy  of  low-order  cumulant  approximations,  namely 
mean  and  variance,  under  cumulant -neglect  moment  closure 
methods  will  approach  the  exact  values  as  the  level  of 
truncation  goes  to  infinity.  Since  the  computational 
effort  is  positively  related  to  the  level  of  truncation, 
however,  our  objective  was  to  find  a  policy  that  minimizes 
the  level  of  truncation  while  maintaining  reasonable 
approximation  accuracy.  Our  mathematical  and  empirical 
findings  suggest  that  there  are  two  such  'logical' 
truncation  levels. 

In  particular,  let  "s"  be  the  highest  order  of  the 
intensity  function,  and  "i"  be  the  number  of  cumulants 
which  we  wish  to  estimate.  It  follows  that  "s+i-1"  would 
be  the  smaller  of  the  two  truncation  levels  and  "2s"  would 
be  the  larger  of  these.  The  lower  truncation  level 
accounts  for  all  cumulants  that  are  in  the  associated 
differential  equations  generated  by  the  system,  and  the 
higher  level  includes  those  cumulants  that  are  one  order 
removed.  Truncating  below  the  lower  level  leads  to 
significant  approximation  error,  while  above  the  upper 
leads  to  unnecessary  computational  effort. 

We  present  numerical  results  to  summarize  the  above 
truncation  scheme  for  High,  Medium  and  Low  traffic 
situations  in  Figures  3-8.  In  all  in  each  case,  the 
maximum  order  of  the  arrival  and  service  intensities  were 
quadratic  and  our  objective  is  to  estimate  the  Mean  and 
Variance  of  the  system.  Hence,  it  follows  that  the  minimum 
truncation  level  is  (2+2-1) =3,  and  the  Maximum  Truncation 
level  is  (2*2) =4.  We  observe  that  mean  is  quite  well 
approximated  under  the  Low,  Medium  and  High  traffic 
intensities.  The  approximation  accuracy  of  the  variance, 
however  varies.  In  particular,  the  variance  is  under¬ 
approximated  under  the  lower  truncation  limit  and  over¬ 
approximated  under  higher  truncation  limit  under  the 
varying  traffic  intensities. 
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High  Traffic 


Figure  3:  Mean  Approximations  under  High  Traffic 
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Figure  4:  Variance  Approximations  under  High  Traffic 
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Figure  5:  Mean  Approximations  under  Medium  Traffic 
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Figure  6:  Variance  Approximations  under  Medium  Traffic 
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Figure  7:  Mean  Approximations  under  Low  Traffic 
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Figure  8:  Variance  Approximations  under  Low  Traffic 
Effect  of  Traffic  Intensities  on  Approximation  Accuracy 

Investigations  into  optimal  truncation  policies  lead  to  the 
observation  that  the  traffic  intensity  greatly  affects  that 
accuracy  of  the  moment  closure  approximations.  In 
particular,  empirical  evidence  indicates  that  the  accuracy 
of  the  approximations  is  best  for  medium  traffic  systems 
and  deteriorates  as  we  move  to  low  and  high  traffic 
systems.  This  finding  is  dependent  on  the  systems  under 


study,  yet  holds  for  a  large  body  of  models.  A  graphically 
illustration  of  this  is  given  in  Figure  9.  To  study  the 
effect  of  traffic  intensity,  we  consider  a  single  node 
system  with  2nd  order  polynomial  intensities.  Generating 
exact  solutions  for  the  low- order  order  cumulants  of  the 
system  via  Kolmogrov  equations,  we  noticed  that  in  general 
the  4th  order  cumualant  was  large  under  low  traffic,  small 
under  medium  traffic,  and  very  large  under  high  traffic. 

It  follows  that  under  3rd  order  truncation,  this  potentially 
large  term  may  be  missing  from  the  equation.  The  effect  of 
this  and  other  large  higher-order  cumulants  may  propagate 
through  as  the  truncation  level  is  increased.  This  aspect 
of  the  project  remains  in  the  exploratory  stages,  yet  is  an 
area  that  warrants  future  research  attention. 


Figure  9:  Traffic  Intensity  and  Approximation  Error 


Analysis  of  the  Sortie  Generation  Process 

Moment  Closure  methods,  under  a  cumulant -neglect  policy, 
were  applied  to  the  sortie  generation  model.  The  model 
considered  mimicked  that  of  Dietz [1] .  A  paper  on  this 
subject  was  written  and  submitted  to  Naval  Research 
Logistics,  yet  is  presently  under  revision  for  a 
resubmission.  The  original  draft  of  this  paper  is  included 
in  Appendix  D,  in  which  a  Phase-type  distribution  is  used 
to  model  the  fork- join  nodes. 

RESEARCH  MATERIALS 

The  following  research  publications  were  developed  under 
this  grant. 


•  Jayaraman,  R.,  Matis,  T.  and  Guardiola,  I.  (2004)  "Effect 
of  Polynomial  Intensity  Functions  on  Cumulant  Derivation 
Procedures",  Proceeding  of  the  2004  Industrial  Engineering 
Research  Conference . 

•  Matis,  T.  I.  and  Kharoufeh  J.  P. ,  (2005)  "Transient 

Queueing  Network  Analysis  of  Sortie  Generation"  Naval 
Research  Logistics  (under  revision) . 

•  Matis,  T.  and  Guardiola  I.  (2005)  "On  the  Stability  of 
Moment  Closure  Methods"  Operations  Research  Letters  (in 
preparation) . 

This  work  was  presented  at  several  research  conferences, 
including 

•  Stability  Region  Identification  for  Non-Linear 
Stochastic  Systems  Using  Moment  Closure  Methods, 
Institute  for  Operations  Research  and  the  Management 
Sciences  2004  Annual  Conference,  Denver,  Colorado, 
October  2004 

•  Jayaraman  R.,  Matis  T.  I.,  and  Guardiola  I.  (2004) 
"Effect  of  Polynomial  Intensity  Functions  on  Cumulant 
Derivation  Procedures",  Proceedings  of  the  Institute 
of  Industrial  Engineering  2004  Annual  Conference. 


CONCLUSION 

The  research  performed  under  this  grant  has  lead  to 
significant  advancements  in  the  theoretical  and  applied 
aspects  of  Moment  Closure  methods.  Notable  theoretical 
advancements  include  documented  studies  on  the  stability, 
accuracy,  and  optimality  of  moment  closure  methods,  and 
applied  advancements  include  efficient  computational 
routines  for  analyzing  large  systems  and  demonstration  of 
applicability  to  the  sortie  generation  process.  A  website 
at  http ; //engr . nmsu . edu/~tmat is  contains  copies  of 
computational  codes,  papers,  and  presentation  produced 
under  this  grant. 
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APPENDIX  A 

MOMENT  CLOSURE  MATHEMATICA©  PACKAGE 

(MCP.m) 


********* 


This  file  is  intended  to  be  loaded  into  the  Mathematics 
kernel  using  the  package  loading  commands  Get  or  Needs. 

*********************************************************** 

************) 

(♦This  Notebook  is  the  keeper  of  all  function  which  are 
necessary  in  order  to  \ 

evaluate  a  give  node  system  using  transient  analysis 

corresponding  to  moment  \ 

neglect*) 

(♦Packages  needed  in  order  to  use  this  program*) 

(♦Wirtten  by  Ivan  Guardiola  &Tim  Matis*) 

<<DiscreteMath 'Combinatorics' 
ccmomcum.m 

(♦This  function  build  the  left  hand  side  of  the  of  the 
partial  differential  \ 
equations.*) 

LHSPDE [Nodes, TRUNC] : * 

Module  [{AA,  AB,  AC,  AD,  AE,  AF,  AG,  AH,  AI ,  AJ,  AK,  AL,  AM,  AN,  AO,  AP,  AQ, 
AR, AS, AT, AU, AV, 

ORLI,MGFList ,  CGFlist,  IW,mgf  1,  Cgf ,  MGF,LHSpde}, 

AA=Table [Table [0, {Nodes}] , {Nodes}] ;AB=MapThread[K, Transpose 
[AA]  ]  ; 

AC=MapThread[m, Transpose [AA] ] ; Clear [AA] ; 

AD=AB/ . {K [b _ ] \ [Rule] K[b] [t] };AE=Take [First [AD] ] ;AF=Take [Fi 

rst  [AC]  ]  ; 

Evaluate [AE] ==0; Evaluate [AF] «l;AG={a,b,  c,  d, e,f,g,h,r,s,t,u 
, v,w,x,y, z} ; 

AH=Take [AG, Nodes] ;AI=Table [0, {Nodes}] ;AJ=Table [TRUNC, {Nodes 

}]; 

AK=Transpose  [{AH,AI,AJ}]  ; AL=MapThread  [m,AH,0]  ;AM=MapThread[ 
K,  AH,  0]  ; 

AN=Table [Subscript [\ [Theta] ,i] , {i, Nodes}] ;AO=ANAAH;AP=Times 
@@A0; 


AQ= Factorial [AH] ; AR=Times@@AQ; 


AS=AP/AR;AT=AM/. {K [z _ ] \ [Rule]K[z] [t] };AU=AL*AS;AV=AT*AS; 

(m[i _ ] / ;Plus@@{i}\ [Great erEqual] TRUNC+1) = 

0; (K[i _ ] [t] / ; Plus@@{i}\ [GreaterEqual] TRDNC+1) =0; 0RLI={} ; 

For [i=l, icTRUNC+l, 

If [i\ [Equal] l,AppendTo [ORLI, Compositions [i, Nodes] ] , 
AppendTo [ORLI , Join [ORLI  [  [i - 
1] ] , Compositions [i, Nodes] ] ] ] ; i++] ; 

order Li s  1 1 =Las  t [ORLI ] ; ORLI = { } ; MGFLi s  t = {  } ; CGF1 is  t = { } ; 

For [i=l,i<Length[orderListl] +1, 

IW=Thread  [AH\  [Rule]  orderListl  [  [i]  ]  ]  ; 

AppendTo  [MGFList,  AU/ .  IW]  ; 

AppendTo  [CGFlist,AV/.IW]  ;i++]  ;mgfl=l+Plus@@MGFList; Cgf=Plu 
s@@CGFlist; 

MGF=Exp [Cgf] ; 

LHSpde=D [MGF, {t, l}] / . {\ [ExponentialE] ACgf 
\  [Rule]mgfl}; {LHSpde,  MGF, Cgf, 
mgf  1} 

] 

(♦This  function  builds  and  evaluates  all  necessary 
symbolics  in  order  to  \ 

create  the  right  hand  side  of  the  partial  differential 
equations*) 

EQUATE [F  , B_, TRUNC] : = 

Module [ {F2 , Y , Y2 , Y5 , Y7 , Nodes , AN, AG, AH, ord, Y8 , Yll , ORLI , orderL 
istl, 

relationsmod, relations, momcumrule, AO, Tl, Y12 , LHL, LHSpde, MGF, 
Cgf ,mgf 1, stu, 

butt,butt2}, F2=Partition [F, 1] ; 

Nodes=Length [B [ [1] ] ] ; 

Y=Table [Join [Take [F2 [ [i] ]] ,Take [B [ [i] ] ] ] , {i, 1, 

Length [F] }] / . {Plus\ [Rule] List} ; 

Y2=Partition [ 

Flatten [Table [ 

If  [ListQ [Y[ [i, 1] ] ] \ [Equal] False, Y [ [i] ] , 

Table [{Flatten [Y [ [i,l] ] J  [ [j] ] ,Drop [Y [ [i] ] , 1] }, { j , 1, 

Length [Y [ [i, 1] ] ] }]  ] , {i, 1, Length [Y] }] ] , 1+Length [B [  [1] ] ] ] ; 


Y5=Y2/ .Thread [Variables [F] \ [Rule] Table [1, {Length [Variables  [ 
F]  ]  }] ]  ; 

AN=Table [Subscript [\ [Theta] , i] , {i, Nodes}] ; 

AG={a,  b,  c,  d,  e,f,g,h,r,s,t,u,v,w,x,y,  z  } ;  LHL=LHSPDE  [Nodes,  TRU 
NC]  ; 

LHSpde=LHL[ [1] ] ;MGF=LHL [ [2]  ]  ; Cgf =LHL [ [3]  ] ;mgfl=LHL[ [4] ] ;AH= 
Take [AG, Nodes]  ; 

Y7=Table[Y5 [ [i , 1] ] *Sum[  (Plus@@(AH*AN) ) Aj/j ! , { j , 1,TRUNC>] /.T 
hread [ 

AH\ [Rule] Take [Y2 [ [i] ] , {2, 2+Length[B [ [1] ] ] - 
1}]] ,{i,l. Length [Y5]}] ; 

ord=Table [Exponent [Y2 [ [i, 1] ] , Variables [F]  ]  ,  (i,  1, Length [Y2] } 

]  ; 


Y8=Table [Prepend [Table [{AN[ [j] ] , ord [ [i,  j]  ]  },  { j , 1, Length [ord 
[[ill] }] , 

MGF]  ,  { i , 1 , Length [Y5 ] } ] ; 

Yll=Plus@@ (Y7* (Table [D@@Y8 [ [i] ] , {i, 1, Length [Y8] }] ) / . {\  [Expo 
nentialE] A 

Cgf  \ [Rule]mgfl}) ;ORLI={}; 

For [ i = 1 , i <TRUNC+ 1 , 

If [i\ [Equal] l,AppendTo [ORLI , Compositions [i, Nodes] ] , 
AppendTo [ORLI , Join [ORLI [ [i - 
1] ] , Compositions [i, Nodes] ] ] ] ;i++] ; 
orderListl=Last [ORLI] ;ORLI={}; 
relationsmods 

MomCumConvert [#, ForMomentQ\ [Rule] "Y" , CenteredQ\ [Rule] "N" , 

MomentSymbol\ [Rule] m, CumulantSymbol\ [Rule]K] &/@orderListl; 

relations=relationsmod/ . {K [r _ ] \ [Rule] K [r] [t] } ; 

momcumrule=relations/.{Equal\ [Rule] Rule} ; AO=AN^AH; 
Tl=Table [Times®@AO/ .Thread [AH->orderListl [ [i] J] ,{i,l. 
Length [orderListl] }] ; 
butt=Table [Thread [AN- 

>Sign [orderListl [ [i] ]] ] , {i, 1, Length [orderListl]  }] ; 

Y12={}; 

butt2= 

Table  [Table [ 

If [Sign [orderListl [ [ j , i] ] ] \ [Equal] 1, a, butt [ [ j , i] ] ] , {i, 1, 


Length [orderListl [ [1] ]  ]  }]  ,  { j , 1, Length [orderListl] }] ; 

For  [i=l, i<Length [butt2] +1, AppendTo [Y12 , DeleteCases [butt2  [  [i 
]  ]  /  a]  ]  ;  i++]  > 

{LHSpde , Y12 , T1 , Yll , momcumrule } 

] 

(*This  function  does  the  moment  matching  operation*) 
EquateCoef ficients [F_, B_, TRUNC] : = 

Module [{LLL, BF, LHSpde, Y12, Tl, Yll, momcumrule} , 

LLL=EQUATE [F, B , TRONC] ;LHSpde=LLL[ [1]  ]  ; Y12=LLL [ [2]  ]  ;T1=LLL[[ 

3]],- 

Y11=LLL [ [4] ] ;momcumrule=LLL [ [5]  ]  ; 

BF=Flatten[ 

Table [Coefficient [LHSpde/. Y12 [ [i] ] ,T1 [ [i] ) ] \ [Equal] 
Coefficient [Yll/ . Y12 [ [i] ] , Tl [ [i] ] ] , {i , 1, 
Length [Tl] }] /.momcumrule] 

] 

(♦This  function  set  initail  conditions  if  any  and  produces 
the  Equations  for  \ 

NDSOlve  thus  this  will  allow  the  building  of  all  necesarry 

function  in  order  \ 

to  solve  the  PDEs  numerically*) 

InitialConditions [F_,B_, TRUNC, Node ID_, Initial] : > 

Modul e [ {ORLI , orderLi s t 1 , BE , 

Nodes , BE1 , BE2 , neweqns , neweqmod , BG , BH , BI , B J , BK , 

BL,BM,BF}, 

Nodes ^Length [B [ [1] ] ] ; 

ORLI={}; 

For [i=l,i<TRUNC+l, 

If [i\  [Equal] 1, AppendTo [ORLI, Compositions [i, Nodes] ] , 
AppendTo [ORLI , Join [ORLI [ [i- 
1] ] , Compositions [i, Nodes] J ] ] ;i++] ; 

orderListl=Last [ORLI] ; Clear [ORLI] ;BE=MapThread [K, Transpose [ 
orderListl]  ]  ; 

BE1=BE/ . {K [y _ ] \ [Rule] K [y]  '  [ t] } ; BE2  =BE/ . {K [y _ ] \ [Rule] K  [y]  [ 

0] \ [Equal] 0}; 

BF=EquateCoef ficients [F,B, TRUNC] ;neweqns=Solve [BF, BE1] ; 
neweqmod=neweqns/ . {Rule\ [Rule] Equal} ; 

If  [Length [NodelD] \ [GreaterEqual] l&&Length [Initial] \ [Greater 
Equal] 1,BG={}; 


For  [j=l, j<Length[Initial] +1, 

For [i=l, i<Length [orderListl] +1, 

If [NodelD [ [j] ] \ [Equal] orderListl [ [i] ] , AppendTo [BG, i] ] ; i++] ; 
j++]  ; 

BH=MapThread [K, Transpose [NodelD] ] ;BI=BH/.{K[n _ ] \ [Rule]K[n] 

[0] };BJ={}; 

For  [i=l,i<Length [Initial] +1, AppendTo [BJ,BI [ [i] ] \ [Equal] Init 
ial [ [i] ] ] ; 

i++] ;BK={}; AppendTo [BK, ReplacePart [BE2,BJ[ [1] ] ,BG [  [1] ] ] ] ; 

If  [Length [Initial] >1, 

For [i=l, icLength [Initial] +1, 

If  [i\[NotEqual]l, 

AppendTo [BK, ReplacePart [BK [ [i- 
1]  ]  ,BJ  [  [i]  ]  ,  BG  [  [i]  ]  ]  ]  ]  ;  i++]  ]  ; 

BL=Last [BK] ;BM=Join [BF, BL] , BM=Join [BF, BE2] ; ] ; 

{BM,BE} 

] 

{♦This  function  obtains  the  numberical  evaluated  functions 
and  return  \ 

interpolating  functions*) 

NumericalSolution [BM_, BE_, TimeAxis_] : = 

Module [{rs}, rs=NDSolve [BM, BE, {t, 0 , TimeAxis} , MaxSteps\  [Rule] 
10000]] 

(♦This  function  calls  all  other  functions  as  well  as  it 
runs  all  plotting  *) 

MomentNeglect  [F_,B_,TRUNC_,NodeID_,  Initial_, TimeAxis_,  PlotF 
uncs_, Range List] : = 

Module [{Sol, INT  ,BM,BE>, 

INT=InitialConditions [F, B, TRDNC, NodelD, Ini tial] ;BM=INT[  [1] ] 
;BE=INT [ [2] ] ; 

Sol=NumericalSolution [BM, BE, TimeAxis] ; 

For [isl, i<Length [PlotFuncs] +1, 

If [Length [RangeList] >=1, 

NotebookWrite [NotebookCreate [] , 

Cell [GraphicsData ["PostScript", 

Display-String  [ 

Plot [Evaluate [{Extract [PlotFuncs, i] }/.Sol] , {t, 0, TimeAxis}, 

PlotRange\ [Rule] Extract [RangeList, i] , 
AxesLabel\ [Rule] {t, Extract [PlotFuncs,!] }, 
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1.  Introduction 

The  Moment  Closure  Program  (MCP)  is  a  Mathmatica®  package  that  estimates 
the  cumulants  (moments)  of  a  non-linear  stochastic  system  over  the  transient  period.  The 
program  does  not  assume  a  particular  parametric  state-distribution  for  the  system  apriori, 
which  results  in  the  simple  truncation,  i.e.  equating  to  zero,  of  all  non-estimated 
cumulants.  This  program  was  originally  conceived  by  Dr.  Timothy  Matis  in  1 998  and 
was  refined  to  its  present  form  by  Ivan  Guardiola  and  Karl  Adams  during  the  2003  and 
2004  academic  years.  The  project  was  sponsored  by  grant  #F49620-03-l-0310  from  the 
Air  Force  Office  of  Scientific  Research,  whose  support  made  this  work  possible. 

2.  Before  You  Start 

Before  using  the  MCP  program,  place  the  file  “MOMCUM.m”  and  “MCP.m”  in 
the  folder  where  the  Mathematica  program  and  kernel  are  located.  This  is  typically  in  the 
directory  C:\Program  FilesWolfram  Research\Mathematica\5.0  for  PC  versions  of 
Mathematica.  The  package  “MOMCUM.m”  was  authored  by  Dr.  Qui  Zheng,  presently 
with  the  Department  of  Bio-Statistics  at  Texas  A&M  University,  and  is  necessary  for 
converting  moments  to  cumulants  as  part  of  the  internal  calculations  of  the  “MCP.m” 
package. 

3.  Entering  a  Stochastic  System 

The  first  step  in  entering  a  stochastic  system  into  Mathematica©  is  to  read  in  the 
MCP  package.  This  is  accomplished  by  typing  the  following  command  into  the  first  line 
of  a  notebook. 


«MCP.m 


On  the  second  line,  a  function  that  calls  information  in  this  package  is  entered.  The 
arguments  of  this  function  describe  all  necessary  information  about  the  network.  This 
function  is  given  below,  and  the  arguments  are  described  in  the  following  subsections. 
MomentNegIect[F,B, TRUNC,NodeII), Initial, TimeAxis,PlotFuncs,RangeList] 

3. 1.  Entering  the  Intensity  Functions 

Instantaneous  changes  in  the  state  of  the  system  in  (t,t+At)  and  their 
corresponding  intensity  functions  are  entered  into  the  lists  B  and  F  respectively.  It 
follows  that  the  lists  B  and  F  must  be  of  the  same  length  and  that  the  elements  of  B  and  F 
correspond  by  position.  Note  that  B  is  entered  as  a  list  of  lists  and  F  is  entered  as  a  list  of 
functions. 

The  variables  of  the  intensity  functions  must  be  entered  as  a  subscripted  letters 
(note:  letters  without  a  subscript  will  not  work  properly).  As  an  example,  the  variable 
corresponding  to  the  state  of  the  first  node  should  be  entered  as  xj.  All  other  elements  of 
the  intensity  function  should  be  entered  in  numeric  form.  The  following  two  examples 
consider  systems  that  are  subject  to  1)  unit  changes  and  2)  bulk  changes  in  (t,t+At). 

3.1.1.  Systems  with  Unit  Changes 

Consider  a  2-node  system  with  instantaneous  unit  changes  in  the  state  of  the 

system.  In  particular,  arrivals  to  the  first  node  occur  at  a  rate  3,  departures  from  the  first 
node  and  the  subsequent  arrivals  to  the  second  node  occur  according  to  the  rate  function 
4  xi2,  and  departures  from  the  second  node  occur  according  to  the  rate  function  5  X23.  To 
specify  this  system,  the  following  lists  would  be  defined  for  B  and  F. 

B={{1,0}, {-1,1}, {0,-1}} 

F={3,4x,2,5x23} 


3.1.2.  Systems  with  Bulk  Changes 

Consider  a  single  node  system  with  bulk  arrivals  and  unit  departures.  In 

particular,  arrivals  to  the  first  node  occur  3  at  a  time  at  a  rate  5  and  departures  occur 
according  to  the  rate  function  4  xi2.  To  specify  this  system,  the  following  lists  would  be 
defined  for  B  and  F. 

F={5, 4  x,2} 

3.2.  Defining  the  Truncation  Level 

The  first  element  defined  by  the  user  is  the  truncation  level  of  the  cumulants. 

This  value  may  be  set  arbitrarily  high  within  the  limits  of  available  memory,  yet  is 
typically  2,3,  or  4  as  this  is  common  for  many  moment  closure  problems.  An  example  of 
the  truncation  level  set  to  3  is  given  below. 

TRUNC=3 

3.3.  Defining  the  Initial  Conditions 

Next,  the  initial  conditions  are  specified.  There  are  three  conditions  that  may 

exist  for  the  system,  1)  all  cumulants  set  to  zero,  2)  one  non-zero  cumulant,  and  3)  more 
than  one  non-zero  cumulant. 

3.3.1.  All  Cumulants  Set  to  Zero 

The  default  initial  condition  of  the  program  sets  all  cumulants  to  zero  at  time  zero. 
Hence,  leave  the  lists  “NodelD”  and  “Initial”  empty  as  shown  below. 

NodeID:={} 


Initial:— {} 


3.3.2.  One  Non-Zero  Cumulant 

If  the  initial  conditions  consist  of  one  non-zero  cumulant,  identify  that  cumulant 
in  “NodelD”  and  enter  the  non-zero  value  in  “Initial”.  Note  that  “NodelD”  is  entered  as 
a  list  of  lists  and  “Initial”  is  entered  as  a  list  of  numbers.  As  an  example,  suppose  at  time 
=  0  there  are  five  items  in  node  three  of  a  five  node  system  and  that  all  other  nodes  are 
known  to  be  empty.  To  enter  this  initial  condition,  specify  “NodelD”  in  the  form  of  a  list 
using  1  in  the  third  position  to  indicate  the  first  cumulant  of  the  third  node  and  then 
specify  5  in  the  first  position  of  the  “Initial”  list  as  shown  below. 

NodeID:={ {0,0, 1,0,0}} 

Initial:={5} 

3.3.3.  Multiple  Non-Zero  Cumulants 

The  procedure  for  entering  initial  conditions  for  multiple  non-zero  cumulants  is  a 

simple  extension  of  that  for  one-cumulant  only.  The  non-zero  cumulants  are  identified  as 
a  list  in  “NodelD”  and  their  values  are  entered  in  the  corresponding  position  in  “Initial”. 
As  an  example,  suppose  at  time  =  0  it  is  known  that  there  are  1 0  items  in  node  three  and 
that  the  expected  number  of  items  in  node  4  is  12  with  a  variance  of  5.  These  initial 
conditions  are  entered  as  shown  below. 

NodeID:={{0,0, 1,0,0},  {0,0, 0,1,0}, {0,0, 0,2,0}} 

Initial:={10,12,5} 

3.4.  Defining  the  Output  Plots 

The  variables  that  control  the  output  plots  of  the  network  are  “PlotFuncs”, 

“TimeAxis”,  and  “RangeList”.  The  cumulants  that  the  user  would  like  to  view  are 
entered  into  the  list  “PlotFuncs”  as  the  elements  K[ij,...][t],  where  i,j, _ are  the  non¬ 

negative  integers  corresponding  to  this  set  of  cumulants.  Note  that  any  desired  marginal 


and  cross  cumulants  up  to  the  order  of  truncation  may  be  specified  in  this  list.  As  a 
reminder,  the  first  order  marginal  cumulants  correspond  directly  to  expectation,  the 
second  to  the  variance,  and  the  third  to  the  skewness.  Likewise  the  second  order  cross 
cumulants  correspond  directly  to  the  covariance  between  the  respective  nodes.  As  an 
example,  consider  a  fully  defined  5-node  stochastic  network  and  assume  our  interest  lies 
only  in  the  expectation  and  variance  of  the  first  node.  The  specification  of  these  output 
plots  would  be  entered  in  the  list  “PlotFuncs”  as  shown  below. 

PIotFuncs={K[2, 0,0, 0,0,0]  [t]  JK[1 ,0,0,0,0,0]  [t]} ; 

The  range  of  the  x  and  y  axis  of  these  plots  are  specified  through  “TimeAxis”  and 
“RangeList”  respectively.  In  particular,  “TimeAxis”  controls  the  range  of  the 
independent  time  variable  on  the  x-axis  and  “RangeList”  specifies  the  range  the 
cumulant  on  the  y-axis.  “TimeAxis”  does  not  have  a  default  value  and  must  be  specified 
by  the  user  for  the  program  to  function  properly.  “TimeAxis”  is  entered  as  a  single  real 
number  that  denotes  the  end  point  of  the  interval  [0,  “TimeAxis”].  “RangeList” 
defaults  to  an  interval  which  is  large  enough  to  show  the  entire  value  of  the  cumulant  on 
the  y-axis.  Often  “RangeList”  is  left  blank  for  the  first  run  for  the  user  to  visual  the 
entire  graph.  “RangeList”,  is  modified  by  entering  embedded  lists  of  the  form  {lower 
endpt  of  y,  upper  endpt  of  y)  that  correspond  by  position  to  those  plots  specified  in 
“PlotFuncs”.  Hence,  the  lists  “PlotFuncs”  and  “RangeList”  must  be  of  the  same 
length.  As  an  example,  consider  the  two  cumulants  specified  in  “PlotFuncs”  in  the 
previous  example.  Suppose  the  user  would  like  to  solve  for  these  cumulants  from  time 
te  [0,50]  and  plot  K[2,0,0,0,0,0][t]  on  the  interval  {0,25}  and  plot  K[l,0,0,0,0,0][t]  on  the 


interval  {5,50}.  The  values  of  “TimeAxis”  and  “RangeList”  would  be  entered  as  shown 
below. 

TimeAxis=50; 

RangeList={{0,25},{5,50}} ; 

4.  Running  the  MCP  Program 

The  MCP  program  may  be  run  by  clicking  on  the  “Kernel”  menu,  then 
“Evaluation”,  then  “Evaluate  Notebook”.  This  will  generate  the  specified  output  graphs 
labeled  with  the  corresponding  cumulant  in  their  own  window.  For  example,  the  graph  of 
the  first  cumulant  of  the  first  node  will  have  the  label  “K[l,0,0,0,0,0][t]”  and  will  open  in 
the  window  “Untitled_l”. 

5.  Limitations 

The  program,  as  written,  has  a  limit  of  twenty  six  nodes.  If  the  system  under 
consideration  contains  more  than  twenty  six  nodes,  the  number  of  variables  in  a  list  that 
is  embedded  is  in  the  internal  calculations  section  of  the  program  must  be  increased.  To 
do  this,  find  the  list  “AG”  in  the  “MCP.m”  pacakage  and  add  variables  to  this  list. 

6.  Contact  Information 

For  further  information  about  the  MCP  program,  please  contact  Dr.  Timothy  I. 
Matis  at  tmatis@nmsu.edu. 


APPENDIX  C 

SAMPLE  IMPLEMENTATION  OF  MCP.m 


«  M3».m 

M=mentlfeglect[{3,  xx  ,  X2},  {{1,  0},  {-1,  1},  {0,  -1}}# 
3,  {},{},  10,  {K[l,  0]  [t] ,  K[0,  l][t]},  {All,  All}] 
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Abstract 

This  paper  develops  a  procedure  for  the  transient  analysis  of  an  aircraft  sortie 
generation  process  that  is  modelled  as  a  closed  network  of  state-dependent  queues. 
The  procedure  uses  the  cumulant  function  method  for  queueing  networks  to  solve  for 
key  time- variant  performance  measures  such  as  the  mean  workload  at  each  station  and 
the  sortie  generation  rate.  Moreover,  the  model  incorporates  a  phase-type  (PH)  service 
time  distribution  for  maintenance  activities  and  accounts  for  aircraft  blocking  at  this 
node.  We  demonstrate  the  implementation  of  the  transient  technique  by  means  of  a 
notional  example. 
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1  Introduction 


In  a  highly  volatile  world  environment,  military  decision  makers  are  currently  faced 
with  the  daunting  task  of  accurately  assessing  the  ability  of  their  forces  to  carry  out  critical 
missions  in  an  efficient  and  expedient  manner.  The  successful  execution  of  such  missions, 
particularly  those  involving  military  aircraft,  hinges  upon  the  availability  of  resources  such  as 
personnel,  aircraft,  munitions,  and  maintenance  facilities.  Assessing  the  operational  readi¬ 
ness  of  a  given  air  base  is  therefore  of  paramount  importance  to  military  decision  makers. 

The  concept  of  operational  readiness  can  refer  to  a  number  of  performance  measures 
within  this  context.  One  important  measure  is  that  of  resource  availability  which  refers  to 
the  proportion  of  time  the  air  base  is  able  to  provide  all  necessary  resources  and  personnel 
to  perform  a  mission.  A  few  examples  of  these  resources  include  diagnostic  equipment, 
tools,  replacement  aircraft  components,  and  physical  hangar  space  necessary  to  perform 
maintenance  activities.  Ultimately,  decision  makers  are  interested  in  the  number  of  successful 
missions  flown  over  a  critical  period  of  time.  This  measure  is  often  referred  to  as  the  aircraft 
sortie  generation  rate.  It  is  obvious  that  the  sortie  generation  rate  is  highly  dependent 
on  the  flow  of  ground  operations.  Perhaps  the  most  critical  portion  of  the  overall  sortie 
generation  process  is  the  maintenance  activity.  This  is  due  to  the  fact  that  distinct  aircraft 
may  have  vastly  differing  maintenance  requirements  upon  sortie  completion.  For  this  reason, 
it  is  crucial  to  assess  such  measures  as  the  current  workload  at  the  maintenance  station(s) 
as  well  as  resource  availability.  However,  these  measures  are  not  necessarily  time  invariant, 
particularly  in  wartime  scenarios  in  which  prevailing  theater-level  dynamics  may  govern 
ground  operations.  The  time-variant  behavior  of  these  measures  ultimately  impacts  the 
ability  of  an  airfield  to  fly  aircraft  sorties. 

Owing  to  the  inherently  complex  interactions  between  resources  required  for  the  sortie 
generation  process,  analysts  have  typically  employed  computer  simulation  techniques  for  the 
purpose  of  evaluating  operational  readiness.  In  the  United  States  Air  Force,  for  instance, 
some  typical  models  that  have  been  employed  are  the  Logistics  Composite  Model  (LCOM) 
[3]  and  the  Sortie  Generation  Model  (SGM)  [1].  Though  such  simulations  allow  analysts 
to  assess  the  utilization  of  resources  (and  other  performance  metrics),  the  implementation 
of  such  models  can  be  cumbersome  due  to  extensive  data  input  requirements  and  long  run 
times  for  a  single  replication.  These  problems  are  exacerbated  in  a  real-time  setting  when 
decision  makers  need  reasonable  answers  in  an  expedient  manner. 
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In  order  to  address  the  shortcomings  of  simulation  models  of  air  base  operations,  some 
authors  have  proposed  analytical  models  to  measure  some  basic  aggregate  features  that  can 
be  used  to  quickly  and  adequately  answer  important  questions  regarding  operational  readi¬ 
ness.  More  specifically,  Dietz  and  Jenkins  [2]  provided  a  mathematical  framework  to  address 
the  problem  of  modelling  the  theater-level  dynamics  of  the  aircraft  sortie  generation  process 
as  a  closed  queueing  network.  In  that  work,  the  authors  presented  a  steady-state,  mean  value 
analysis  (MVA)  for  several  important  performance  measures  involved  in  the  aircraft  sortie 
generation  process.  The  key  innovation  in  their  model  was  the  incorporation  of  a  single, 
fork-join  node  in  the  queueing  network  that  enables  the  analytical  modelling  of  concurrent 
maintenance  activities  subsequent  to  sortie  completion.  This  approach  to  the  problem  is 
significant  for  several  reasons.  First,  it  allows  for  an  aggregation  of  many  complexities  into 
single-  or  multi-server  queueing  stations.  Second,  it  provides  an  analytical  framework  upon 
which  to  build  models  of  higher  resolution  if  needed,  and  third,  it  provides  fast  numerical 
results  for  decision  makers  who  require  a  “snapshot”  of  their  current  operational  capability 
by  avoiding  time-consuming  simulation  replications.  Hackman  and  Dietz  [6]  extended  the 
preliminary  work  of  [2]  by  allowing  the  service  times  at  each  node  of  the  network  (i.e.,  at  each 
stage  of  the  sortie  generation  process)  to  be  arbitrarily  distributed  rather  than  exponentially 
distributed. 

Though  these  two  papers  provide  a  framework  upon  which  the  sortie  generation 
process  may  be  analyzed,  the  work  does  suffer  an  important  shortcoming.  Both  works  provide 
a  steady-state,  mean  value  analysis  (MVA)  of  the  workload  at  each  node  and  the  sortie 
generation  rate  as  opposed  to  a  more  realistic  transient  analysis.  Although  the  steady-state 
analysis  allows  for  mathematical  tractability  of  the  model,  it  fails  to  incorporate  the  realistic, 
time-variant  behavior  of  the  important  measures  of  operational  readiness.  Additionally, 
the  queueing  network  models  developed  in  these  papers  do  not  account  for  the  blocking 
of  aircraft  that  occurs  at  the  maintenance  node  due  to  limited  hangar  space.  Hence,  the 
primary  objective  of  our  work  is  to  extend  that  of  Dietz  and  Jenkins  [2]  and  Hackman  and 
Dietz  [6]  by  explicitly  incorporating  the  time  dependence  of  queueing  performance  measures 
for  a  more  accurate  assessment  of  operational  readiness  by  using  a  new  analytical  technique 
for  the  transient  analysis  of  queueing  networks. 

The  main  contributions  of  this  work  can  be  summarized  as  follows.  We  provide  a 
transient  analysis  for  a  closed  queueing  network  model  of  the  sortie  generation  process.  The 
primary  measure  of  operational  readiness  is  the  time-variant  sortie  generation  rate,  though 
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we  also  consider  the  workload  (or  congestion  level)  of  the  system  measured  by  the  expected 
number  of  aircraft  present  at  each  station  at  a  given  point  in  time.  To  that  end,  we  formulate 
a  phase-type  approximation  for  the  distribution  of  aircraft  holding  time  at  the  fork-join  repair 
node,  and  adapt  and  employ  the  cumulant  function  method  previously  applied  to  queueing 
networks  by  Matis  and  Feldman  [9].  Our  approach  to  the  problem  allows  us  to  compute  the 
first  moment  of  the  aforementioned  measures  as  explicit  functions  of  time. 

The  remainder  of  the  paper  is  organized  in  the  following  manner.  Section  2  reviews 
the  important  features  of  the  steady-state  model  of  [2].  In  section  3,  we  present  our  modified 
queueing  network  model  of  the  aircraft  sortie  generation  process.  Section  4  discusses  the 
cumulant-based,  transient  analysis  of  the  queueing  network  while  section  5  provides  a  nu¬ 
merical  example  demonstrating  the  implementation  of  the  procedure.  Finally,  we  give  some 
concluding  remarks  in  section  6. 

2  Queueing  Network  Model 

Dietz  and  Jenkins  [2]  were  apparently  the  first  to  present  a  formal  mathematical  model 
of  the  sortie  generation  process.  Before  describing  our  transient  analysis  and  extensions  of 
their  model,  we  provide  a  brief  review  of  the  latter.  There  are  six  activities  that  can  be 
identified  in  the  process  of  flying  sorties  from  a  given  air  base.  In  the  queueing  network 
framework  of  [2],  each  activity  is  modelled  as  an  individual  queueing  station.  The  flow 
entities  of  the  network  are  a  fixed  number  of  aircraft  ( N )  that  pass  through  the  six  nodes 
(stations)  according  to  a  stochastic  routing  matrix  P  :=  [p,j]  where  py  denotes  the  stationary 
probability  that  an  aircraft  completing  activity  i  proceeds  next  to  activity  j.  The  network  is 
assumed  to  be  closed  in  that  the  N  aircraft  never  leave  the  system,  nor  do  additional  aircraft 
arrive  to  the  system.  The  closed  queueing  network  model  consists  of  the  following  nodes: 
pre-flight,  sortie,  troubleshoot,  a  fork-join  maintenance  node,  turnaround,  and  munitions 
upload.  The  fork-join  node  consists  of  five  substations  representing  five  critical  systems  of 
an  aircraft  that  may  or  may  not  require  maintenance  subsequent  to  sortie  completion  and 
the  troubleshooting  activity.  Each  station  is  assumed  to  have  an  exponentially  distributed 
service  time.  Figure  1  gives  a  graphical  depiction  of  the  flow  of  aircraft  through  this  network. 

The  closed  queueing  network  model  was  analyzed  using  a  mean  value  analysis  (MVA) 
heuristic  developed  by  Rao  and  Suri  [11]  to  accommodate  the  nuance  of  a  fork-join  node  and 
enhanced  by  Jenkiris  [7]  to  handle  fork-join  nodes  with  probabilistic  branching.  While  the 
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Figure  1:  Queueing  network  model  of  sortie  generation  process. 


MVA  produces  exact  results  for  product-form  networks,  the  heuristic  MVA  procedure  pro¬ 
vides  only  approximate  performance  measure  values  for  any  closed  network  with  more  than 
one  customer  in  the  system.  Using  this  approach,  the  primary  (steady-state)  performance 
measures  computed  in  [2]  were:  i)  the  mean  number  of  aircraft  at  each  station  (workload), 
ii)  the  throughput  at  each  station,  and  iii)  the  overall  throughput  of  the  network,  which 
directly  corresponds  to  the  time-invariant  sortie  generation  rate.  These  results  are  useful 
due  to  their  mathematical  tractability  and  ease  of  implementation;  however,  they  do  not 
account  for  the  inherently  time- variant  behavior  of  the  station  (or  system)  workload  and  the 
sortie  generation  rate.  Moreover,  the  model  of  [2]  ignores  the  effects  of  blocking  by  assuming 
infinite  capacity  queues  at  each  station;  however,  limited  hangar  space  is  a  reality  at  most 
air  bases  causing  aircraft  to  be  blocked  at  the  troubleshoot  node  of  the  network. 

In  the  following  section,  we  present  a  modified  version  of  the  basic  model  of  Dietz 
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and  Jenkins  [2]  in  order  to  consider  the  transient  (time-dependent)  behavior  of  the  expected 
workload  in  the  system  and  the  sortie  generation  rate.  In  the  subsequent  section,  we  describe 
the  formal  procedure  for  analyzing  the  model  in  the  transient  regime. 

3  Modified  Queueing  Network  Model 

The  first  significant  difference  in  our  model  is  the  characterization  of  the  holding 
time  in  the  maintenance  hangar  as  a  phase-type  (PH)  distribution.  PH  distributions  are 
extremely  useful  when  modeling  the  distributions  of  nonnegative  random  variables.  Figure 
2  gives  a  graphical  depiction  of  our  modified  queueing  network  model. 


Figure  2:  Sortie  generation  model  with  PH-distributed  repair  times. 

It  is  important  to  note  the  distinction  between  the  models  of  Figures  1  and  2.  In  Figure  2, 
it  is  compulsory  to  add  node  6  to  the  maintenance  node  (node  5  of  Figure  1)  in  order  to 
implement  our  two-phase  holding  time  distribution. 
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For  our  modified  model,  we  adopt  the  notation  in  [2].  Referring  to  Figure  1,  let  s,- 
denote  the  mean  service  time  at  node  i,  and  let  r,  denote  the  number  of  servers  at  node  i,  for 
i  =  1,2, 3, 4.  Owing  to  the  fact  that  service  rates  are  state  dependent,  whenever  n  customers 
are  present  at  node  i,  the  state-dependent  service  rate,  //»(«),  is  given  by 

Ati(n)  =  min{n/sj,  n/sj,  i  =  1, 2, 3, 4.  (1) 

With  regard  to  node  5  of  Figure  1,  the  fork-join  repair  node,  we  define  s5>fc  as  the  mean 
service  time  at  repair  station  k,  r3tk  denotes  the  number  of  servers  at  repair  station  k,  and 
g5)fc  is  the  probability  that  repair  activity  k  is  required  by  an  arriving  aircraft,  k  =  1, 2, 3, 4, 5. 
The  state-dependent  service  rates  for  the  maintenance  activities  are  thus 

/z5,fc(n)  =  min{n/s5,jfc,  r5)fc/s5ifc},  k  =  1, 2, 3, 4, 5.  (2) 

The  service  rates  of  Equation  (2)  correspond  to  the  rate  parameters  of  the  associated  expo¬ 
nential  distributions  at  each  of  the  five  repair  stations  of  node  5.  For  our  particular  queueing 
network  model,  the  parameter  values  are  summarized  in  Table  1. 


Table  1:  Model  parameter  values  for  sortie  generation  queueing  network. 


Node  Index 

Mean  Service  Time 

Repair  Probability 

Number  of  Servers 

1 

Si  =  0.25 

7*!  =  OO 

2 

s2  =  2.00 

r2  =  oo 

3 

s3  =  1.25 

r3  =  oo 

4 

s4  =  0.50 

r4  =  1 

5,1 

Ssj  =  2.20 

9s, i  =  0.17 

r5,i  =  1 

5,2 

$5,2  =  2.27 

<75,2  =  0.39 

<N 

II 

<N 

ua 

5,3 

$5,3  =  2.37 

<75,3  =  0.21 

r5,3  =  2 

5,4 

S5>4  =  1.50 

95,4  =  0.27 

T— < 

II 

£ 

5,5 

«5,1  =  1.19 

95,5  =  0.46 

r5, 5  =  2 

It  should  be  noted  that  the  activity  labelled  “munitions  upload”  in  [2]  is  contained  in  the 
“turnaround”  activity  in  this  model  so  that  the  total  mean  service  time  at  the  turnaround 
queue  is  1.25  time  units. 

We  examine  the  operational  readiness  of  a  single  air  base  that  services  N  aircraft  with 
a  maintenance  hangar  capacity  of  2  planes  (i.e.,  we  relax  the  assumption  that  the  mainte¬ 
nance  hangar  has  infinite  capacity).  An  airplane  undergoing  service  will  occupy  a  hangar 
space  until  all  required  service  activities  have  been  completed.  The  troubleshooting  activi¬ 
ties  of  node  4  will  be  halted  when  the  hangar  is  full  and  will  resume  only  when  a  departing 
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aircraft  frees  capacity  in  the  subsequent  repair  node.  The  repair  activities  of  node  5  are  per¬ 
formed  concurrently  and  assumed  to  be  mutually  independent.  In  the  following  subsection, 
we  describe  the  means  by  which  the  original  fork-join  construction  of  the  repair  node  is  used 
to  provide  a  phase-type  representation  of  the  maintenance  holding  time  distribution.  The 
parameters  of  the  phase-type  distribution  will  ultimately  be  used  to  describe  the  stochastic 
evolution  of  the  system,  thereby  allowing  for  our  transient  analysis  procedure. 

3.1  PH  Distribution  for  Repair  Holding  Times 

The  distribution  of  aircraft  holding  time  at  the  fork-join  repair  node  shall  be  repre¬ 
sented  by  a  continuous  three-parameter,  two-phase,  state-dependent  distribution.  We  first 
briefly  review  PH  distributions  which  are  described  in  detail  in  the  pioneering  book  by  Neuts 
[10].  Consider  a  Markov  process  on  a  finite  state  space  £  :=  {1, 2, . . . ,  u-f  1}  such  the  state  j 
is  transient  for  1  <  j  <  u,  and  state  u  + 1  is  the  absorbing  state.  The  infinitesimal  generator 
matrix  Q  of  this  process  may  be  partitioned  as 


where  the  elements  of  the  u  x  u  matrix  G  are  such  that  Gu  <  0  and  Gy  >  0,  i  ^  j.  Moreover, 

Ge  +  G°  =  0 

where  e  is  a  column  vector  of  ones.  The  initial  distribution  of  the  Markov  process  is  given 
by  the  row  vector  a'  :=  (a,  ctu+1)  where  a  is  a  row  vector,  au+i  is  the  probability  that  the 
process  begins  in  state  u  +  1,  and 

ae  +  ou+i  =  1. 

A  probability  distribution  F  for  a  nonnegative  random  variable  is  said  to  be  a  PH  distribution 
if  and  only  if  it  is  the  distribution  of  the  random  time  until  absorption  of  the  aforementioned 
Markov  process.  In  such  case,  the  distribution  is  said  to  have  representation  (a,G). 

The  representation  of  repair  holding  times  necessarily  incorporates  the  state  depen¬ 
dence  of  service  rates  for  individual  repair  activities.  Following  the  notational  convention 
of  Neuts  [10],  the  PH  distribution  is  denoted  by  the  ordered  pair  (a,G(n))  where  G(n)  is 
the  matrix  of  state-dependent  transition  rates.  In  this  work,  the  holding  time  distribution 
at  the  repair  node  has  the  representation 

(S,G(n))  =((1,0).  [  -(-(")  +  ■*<»))  (4) 
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where  i/*  (n) ,  i  =  1,2,3,  denote  the  unknown  parameters  of  the  phase-type  distribution. 
In  order  to  describe  the  stochastic  evolution  of  this  system,  our  aim  is  to  determine  these 
unknown  parameters  by  matching  moments. 

Referring  to  Figure  1,  let  T5<k,  k  G  {1, . . . ,  5},  be  an  exponentially  distributed  random 
variable  representing  the  aircraft  holding  time  at  repair  node  k,  and  let  N5(t)  G  {0,1,2} 
denote  the  number  of  aircraft  in  the  repair  hangar  of  node  5  at  time  t.  The  moments 
E[T£h  |  N5(t)  =  n],  j  =  1, 2, 3,  may  be  calculated  for  repair  activity  k  using  the  exponential 
assumption  of  individual  repair  times  by 

r°° 

E[T’k  \  Nh[t)  =  n]= 

Jo 

Let  Q  —  {0,  {1}, . . . ,  {1, 2, 3, 4, 5}}  be  the  set  of  all  possible  repair  activity  combinations 
for  an  aircraft.  From  the  independence  of  the  repair  activities,  the  probability  that  the  set 
u  G  D  is  required  by  an  arriving  aircraft,  is 

= n  n*1- qik (5) 

Denote  by  a  continuous  random  variable,  T5,  the  total  time  spent  by  an  aircraft  in  the  hangar 
at  node  5  until  all  required  repair  activities  have  been  completed.  The  conditional  moments 
E\Tl  |  N5(t)  —  n],  j  =  1, 2, 3,  may  be  calculated  as 

E[T£  |  N5(t)  =  n]  =  Y^  *uE[raax.  {T£k}  |  N5(t)  =  n].  (6) 

wen 

For  the  queueing  network  model  defined  in  Table  1,  the  moments  of  T5  are  given  by 


Table  2:  Repair  time  moments. 


n 

E[T5  |  N5(t)  =  n] 

E\Ti  |  Ns(t)  =  n] 

E{TH  |  7m  =  "] 

0 

0 

0 

0 

1 

2.07944 

9.33339 

61.6057 

2 

1.82659 

7.63503 

48.3688 

Suppose  X  is  a  nonnegative  random  variable  with  PH  representation  (a,  G).  Then  it  is  well 
known  (cf.  Kao  [4])  that  the  jth  moment  of  X  can  be  obtained  by 

E[X*)  =  (-l)^bG^e,  j  >  1.  (7) 
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For  j  =  1,2,3  and  each  n  E  {0,1,2},  the  moments  of  the  phase-type  distribution 
defined  in  Equation  (4)  are  calculated  and  matched  to  those  given  in  Table  2  using  Eq.  (7) 

by 


^|7V5(t)  =  n]  =  (-l)^!(l,0) 


-(iq(n)  +  i/2(n)) 
0 


-■an;)  « 

This  operation  generates  three  sets  of  three  algebraic  equations  whose  simultaneous  solution 
yields  estimates  for  the  unknown  parameters  i/j(n),  i  =  1, 2, 3,  of  the  phase-type  distribution. 
For  our  defined  queueing  network  model,  these  parameter  estimates  axe  listed  in  Table  3. 


Table  3:  Parameters  of  the  PH  distribution. 


n 

i/i  (n) 

V2  (n) 

"3  (n) 

0 

0 

0 

0 

1 

12.2582 

1.4874 

0.42102 

2 

5.1942 

1.36185 

0.47136 

The  phase-type  representation  of  the  fork-join  nodes  is  incorporated  into  the  model  of  sortie 
generation  depicted  in  Figure  1  yielding  the  model  representation  displayed  in  Figure  2. 
In  what  follows,  we  discuss  the  stochastic  evolution  of  this  closed  queueing  network  before 
presenting  the  transient  analysis  procedure  of  section  4 


3.2  Stochastic  Evolution  of  the  System 

The  evolution  of  this  stochastic,  time- variant  system  is  described  as  follows.  Let  Xi(t) 
denote  the  number  of  aircraft  at  station  i  at  time  t  >  0  and  let  the  multivariate  stochastic 
process  {Xi(t)  :  t  >  0,  i  —  1, 2, . . . ,  6}  describe  the  state  of  the  network  at  time  t.  The  state 
space  Ei  of  Xx(t)  is  given  by  E{  =  {0, 1, 2, ... ,  TV},  for  1  <  i  <  6.  We  will  henceforth  denote 
this  vector- valued,  Markov  process  by  {X(f)  :  £  >  0}  where  X(t)  =  (Xi (t), . . .  ,X6(t)).  In  a 
small  interval  of  time  (t,t  +  6t),  we  note  that,  for  each  node  of  the  queueing  network,  the 
change  of  state  may  decrease  by  one  aircraft,  increase  by  one  aircraft,  or  remain  the  same. 
We  describe  this  collection  of  all  possible  state  changes  by  the  set 


B  =  {(-1, 1, 0, 0, 0, 0),  (-1, 0, 0, 1, 0, 0),  (0,  -1, 1, 0, 0, 0),  (0,  -1, 0, 1, 0, 0),  (0, 0,  -1, 0, 0, 1), 

(0, 0, 0,  -1, 1, 0),  (0, 0, 0, 0,  -1, 1),  (1, 0, 0, 0,  -1, 0),  (1, 0, 0, 0, 0,  -1)} .  (9) 

The  intensity  functions  corresponding  to  the  elements  of  this  set  will  be  denoted  as 

f{b i (■^'1  j  •  •  •  i  "^6) 
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for  Xi  e  Ei.  Combining  information  from  Tables  1  and  3,  the  functions  corresponding  to  (9) 
are  specified  in  Table  4  where  /[*,„]  (z)  is  a  0,1  indicator  function  on  Et  for  x  €  [/,«].  These 
defined  intensity  functions  will  be  represented  as  finite  polynomial  functions  of  the  form 

•  •  •  )  ®6)  0(^1)*'*  >  bfi'jhbiy—  ,be(xi)  >%6)i  (IQ) 

where  <f>{b\,  ■  ■  ■  ,b6)  is  a  constant  and  hbly...yb6(x i,---  ,xe)  is  a  polynomial  function,  for  the 
purpose  of  subsequent  cumulant-based  mathematical  operations. 


Table  4:  Intensity  functions  for  the  aircraft  sortie  model. 


(fell ...ybe)  €  B 

fbi,...ybe  (^1»  *  ♦  •  y%6) 

(-1,1, 0,0, 0,0) 

Pl2(xi/Sl)  =  4.0pi2Xi 

(-1,0, 0,1, 0,0) 

pu{xx/si)  =  i.OpuXi 

(0,-1, 1,0, 0,0) 

P23(X2/S2)  =  0.25P23Z2 

(0,-1, 0,1, 0,0) 

P2i{x2/S2)  =  0.25p24Z2 

(1,0, -1,0, 0,0) 

X3/S3  =  0.802:3 

(0,0, 0,-1, 1,0) 

(2:4/s4)/[0ii](x5  +  x6)  «  2x4  (l  +  0.5  (x5  +  x6)  -  0.5  (x5  +  x6)2) 

(0,0, 0,0, -1,1) 

^(Ol^lo.ojfe)  +  Vl(l)/(l,lj(*5)I’[0,0](a;6)  +  ^1(2)  (/[2,2](*5)f(0,0](a:6)  +  ^[1,1] (^5)^[1,1] (^e)) 

=  21.9193x5  -  9.66108xg  -  7.06402x5x6 

(1, 0,0,0, -1,0) 

^2 (0)/[o, 01(2:5)  +  l^(l)/[l,l](a:5)/|o,o](X6)  +  ^2(2)  (/[2,2)(r5)/[0,0)(r6)  +  /[1,1](X5)/(1,1)(X6)) 

=  2.29388X5  -  0.806475x5  -  .12555xsx6 

(1,0, 0,0, 0,-1) 

•^(O)^[0,0l(x6)  +  ^(Hfll.lJ^eKjO.oK^s)  +  ^3(2)  (/[2,2](2;6)^0,0](a:5)  +  /[l,l](^6)f[l,l](^5)) 

=  0.605482x6  -  0.184451x6  +  .05034x6x5 

The  intensity  functions  of  Table  4  depend  explicitly  on  the  parameters  of  the  PH  distribution 
(Ui(n),  1  <  i  <  3).  In  what  follows,  we  present  the  means  by  which  to  apply  cumulant 
functions  for  the  purpose  of  performing  a  transient  analysis  of  the  Markov  process  {X(<)  : 
t  >  0}.  In  particular,  our  interest  is  in  approximating  the  time- variant  functions  E[Xi{t)\ 
for  i  —  1, . . . ,  6  as  well  as  the  time- variant  sortie  generation  rate. 

4  Cumulant-Based  Transient  Analysis 

There  are  relatively  few  analytical  procedures  through  which  the  Markov  process 
{XM  :  *  >  0}  may  be  practically  analyzed  over  the  transient  period.  Many  of  the  existing 
analytical  techniques  suffer  from  either  a  high  degree  of  mathematical  complexity  or  from 
computational  intractability.  As  an  example,  consider  the  common  approach  of  using  Kol¬ 
mogorov  differential  equations  to  obtain  the  transient  state  probabilities  of  {X(f)  :  t  >  0}. 
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The  number  of  such  equations  is  given  by 


1  +  J2Ci 


(11) 


\  *=1  / 

where  n  denotes  the  number  of  nodes  in  the  network  and  c,-  is  the  capacity  of  the  ith  node. 
As  an  example,  the  previously  defined  queueing  network  with  n  =  6  and  c*  =  N  for  all 
i,  generates  approximately  5.15  x  1010  differential  equations  when  N  =  10  aircraft  and 
7.43  x  1014  equations  when  N  =  50  aircraft.  Hence,  the  number  of  Kolmogorov  differential 
equations  tends  to  intractability  for  even  small  N. 

Alternatively,  a  cumulant-based  approach  may  be  used  to  approximate  the  first  order 
cumulants  (means)  of  the  multivariate  state  distribution  of  (X(t)  :  t  >  0}  in  a  much  more 
computationally  efficient  manner.  Specifically,  the  number  of  differential  equations  using 
these  procedures  is 


(i2) 

where  m  is  the  level  of  cumulant  truncation.  The  specification  of  m  =  2  assumes  that  the 
state  distribution  is  multivariate  normal,  which  is  not  practical  in  most  cases.  It  has  been 
previously  shown,  however,  that  raising  the  truncation  level  to  m  =  3,  which  introduces  a 
skewness  measure,  is  sufficient  for  reasonable  approximations  of  the  first  cumulant  of  Xi{t). 
For  the  defined  queueing  network  {X(t)  :  t  >  0}  with  n  =  6  and  m  —  3,  83  differential 
equations  will  be  generated  by  a  cumulant-based  approach  independent  of  N.  This  set  of 
equations  is  sufficiently  small  to  facilitate  an  efficient  numerical  solution. 

The  use  of  cumulant-based  analysis  procedures  for  a  general  n-node  network  will  first 
be  discussed,  following  which  the  procedure  will  be  applied  to  the  defined  closed  queueing 
network  (X(f)  :  t  >  0}  representing  the  aircraft  sortie  generation  process. 


4.1  General  Procedural  Description 


For  0n,t  >  0  and  1+  the  set  of  non-negative  integers,  let  M[9 1, ..., 6n,t)  be  a 

multivariate  moment  generating  function  for  X(f)  defined  as 

mai.-~.rn1— 9? 


M(e1,...,en,t)  =  T 


(13) 

'a1,...,a„€l+  Oi!---On!  ' 

where  maii...ia„(t)  are  the  joint  moments  of  X(f).  Likewise,  let  K(9ly ...,  9n,  t )  be  a  multivariate 
cumulant  generating  function  defined  as 

=  . -w-6? 

*  n 


Ja\  ,...,a„€l+ 


eql-'-aJ 


(14) 
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where  kait,^an(t)  are  the  joint  cumulants  of  X(f).  By  definition  (Kendall  [5]),  the  generating 
functions  of  Eqs.  (13)  and  (14)  have  the  functional  relationship 


M(9u...,0n,t)  =  eK{dl’’6n't), 


(15) 


through  which  the  joint  moments  define  the  cumulants  of  X(t). 

The  multivariate  moment  generating  function  M{9 1,  ...,0n,t)  and  the  polynomial  in¬ 
tensity  functions  fb1,...,bn(xi,  ...,xn)  of  a  Markovian  network  may  be  related  through  the 
partial  differential  equation 


dM(91,...,9n,t) 

dt 


=  £  (eMi+-+W„_i)/lii  Ai( 

bl,...,bn 


d 

ddi 


bi 


OOi 


6n 


-)M(9u...,9n,t)  (16) 


where  /i>i  (gJ^->  •  •  •  ,  g^-)M(9i, 9n,t)  is  a  partial  differential  operator  that  replaces 
terms  of  the  form  a;’  in  the  polynomial  function  /  with  the  iih  partial  derivative  of  M(9i, ...,  9n,  t) 
with  respect  to  9j.  The  partial  differential  equation  of  Eq.  (16)  may  be  solved  directly  if 
possible,  yet  this  is  often  not  plausible  due  to  the  frequent  complexity  and  non-linearity  of 
the  equation. 

As  an  alternative,  an  mth  order  truncated  cumulant  generating  function,  defined  as 


K°(91,...,9n,t) 


E 


m 


Ol 


0an 


ai,...,a„€Z(m) 


a:! 


(17) 


n 

for  Z(m)  =  {(aj, . . . , an)  :  €  X+,  J2ai  —  m}>  may  be  used  to  extract  a  closed  set  of 

»=i 

ordinary  differential  equations  from  this  partial  differential  equation.  In  particular,  Eq.  (17) 
is  substituted  into  Eq.  (16)  yielding  the  expression 


QeK°(Ou...,en,t) 

dt 


=  Y,  (eb'9'+  -+K9"-l)fbl . &„( 

bi,...,bn 


d 


d 


d9bl  d9br 


> 


K°{e  i— ,0»,t) 


(18) 


Expanding  all  partial  derivatives,  performing  Taylor  series  expansions,  and  substituting  ex¬ 
panded  moment  generating  function  in  Eq.  (18)  generates  the  expression 


E 


(*)0 11  •  •  •  6? 


oi,—  ,On€l+ 

/  oo 


Oi 


!•••«.! 


E 


Ol,— ,a„ez(m) 


y'  |  v~^  (^i^i  + - h  bn9n)k  \  / ^ 

6iti  \h  kl  )' 


“i!  —  on! 

mait...,gn{t)911  ■  -  ■  9“n  \ 
aJ-’-On!  / 


j(bi,...,bn) 


. 6n(^l.-.^n)  , 


(19) 


St 
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where  c6li...,6n(z)  is  a  constant,  j(bu-  •  •  ,6n)  G  1+,  and  Tpllt...tbrt(0i, -,^n)  is  a  polynomial 
function  of  0n  comprised  of  moments  and  non-truncated  cumulants.  Since  both  sides  of 
Eq.  (19)  are  polynomial  expressions  of  0\, ...,  0n,  the  coefficients  of  each  unique  combination 

of  (tfj1  •  •  •  :  Cj  G  l+,ai  -i - (-  an  <  m}  on  the  left  and  right  hand  side  must  equate. 

Hence,  equating  these  coefficients  and  expressing  the  joint  moments  m0l  ,...|0n(<)  as  functions 
of  the  non-truncated  joint  cumulants,  kau...>an(t),  for  {a1(  ■  •  •  an}  G  Z(m)  forms  a  closed  set  of 
ordinary  differential  equations  whose  solution  approximates  the  time- varying,  non-truncated 
cumulants  of  (X(f)  :  t  >  0}  up  to  the  mth  order. 

4.2  Analysis  of  the  Sortie  Generation  Model 

In  this  subsection,  we  apply  the  general  cumulant-based  procedure  to  our  modified 
queueing  network  model  for  the  aircraft  sortie  generation  process.  Let  K°(0 1, ...  ,06,t)  be 
a  third-order  (m  —  3),  truncated  cumulant  generating  function  defined  as  in  Eq.  (17)  for 
the  sortie  queueing  network  model  (X(t)  :  t  >  0}  defined  in  section  3.  Using  Eq.  (18),  the 
function  K°(0i,  ...,06,t)  is  related  to  the  intensity  functions  defined  in  Table  4  through  the 
partial  differential  equation 


o  «■(»„. ..e,,t)  QeK[e ,....06, t)  £  *«>,,..  06,t) 

=  4.0 p12{ee'-9'  -  l)^ss - +  4.0 p14(ee<-0'  -  1) 


8t 


891  '  d9 1 

a  .  deK(9u-eB,t) 

+  0.50p?3(ee3-e>  -  1) - ^ - +0.50 pn(e0*~03  -  1)  ^ - 

a  „  ae/r(e„...9e,t)  a  ir(e1>..-e6,t) 

+  0.80(e®>'®3  -  1)- - — - +  2.0(e®5~®4  -  1)  -gg- - 


( 


x  1  +  0.50 


f  QeK(0i,...6t.i)  deK{0u...ee,t)\  / #2eK(e1,...e<s,t )  d2eK(0t,...0e,t)  02eK(f,,...0„t) 

V  d9l  +  W6  )  “°-50V  '  ^ 


80s  89 5 


+  ■ 


+  (e 


®8-®5  - 1)  (21. 


919388 


QzK[0u  ...ee,i) 

~805 


-  9.66108 


801 

a2eK(9i,...es,t) 

~89\ 


891 


-  7.06402 


+  2°'  805865 

Q2  eK(0l,...0,,t) 


^  // 


(e®3~®5  -  1)  ^2. 

+  (e®3-®6  -  1)  ^0. 


g  K(0u...Os,t)  a2  K(0u...09,t) 

29388 - ^ - 0.806475 - — - 0.12555 


89 586 6 

#2eK(91,...0e,t )> 


89 5 


891 


806  ) 


OeK(01,...0s,t)  a2eK(0u...0s,t) 

605482 - — - 0.184451 - - +  0.05034 


865865 

$2  eK(0i,...0a,t) 


865 


891 


86586, 


St-)’  (20) 


Patterned  after  Eq.  (19),  the  expansion  and  substitution  operations  of  the  previous  subsec¬ 
tion  are  performed  on  Eq.  (20)  yielding  the  general  polynomial  expression 
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( V"  ma i,-.,o6  C0\  f  V''  ^<H o«(*A  _ 

\v^-/0t€X+  ail-’-ae!  /  l  ^~u<h€Z(3)  ai!*»*a6!  J 

4.0pi2  (jr  --  2  ^  ^  ^  ^(-1  ,i,o, 0,0,0) (#i ' '  •  &e)  +  4.0pi4  -  4  1 — ^  ^(-1,0,0, i,o,o)(^i  * •  *^s) 

+  0.50p23  (V^  -  3  ^  ® (0, -1,1,0, o,o)(^i  ’  ‘  ‘  ^e)  +  0.50p24  ^  4— y  ^  1®,(o,-i,o,i,o,o)(®i  •  •  ■< 

+  0.80  tjL-ll)  $(ito_i,o,0,o)(^i -”0e) 


+  0.50p23 


(02  -  0i y 


(h -M 


(#4  -  02) 


-  1  ^(0,-1, 0,1, 0,0' 


I)(#l  -06) 


'®(1,0,-1,0,0,0)(^1  06) 


x  (1  +  0.50  (*go, 0,-1, 1,0) (*1  •  •  •  As))  -  0.50  (*g0i0,_lflj0)(*i  •  *  * *))) 

+  (g<^) 

x  (21.9193^g;o  OiOt_1)1)(0i  •••«»)-  9.661084-g;o  00  •••*)-  7.06402*g>oo  0 >_1  ;1)(0i  •  •  -  06>) 

x  (2.29388*g>,ili0|_li0)(fli  ■■■06)~  O.8O6475*golo_lo)(0,  •  •  •  06)  -  O.12555*g>01>0  _1  o)(0!  •  •  - 06j) 
x  (o.605482®<g0ili0i0i_1)(fl,  --06)-  0.184451^g>0  >1>0  0>_1}(^  •  •  •  06)  +  0.05034*«flfOi0f.1)(91  •  •  •  06))  . 


The  coefficients  of  each  element  of  the  set  {B^B^2  •  •  •  B£6  :  a*  €  T+,  ai  H - h  a6  <  3}  on  the 

left  and  right  hand  side  of  this  expression  are  equated  to  form  a  set  of  83  ordinary  differential 

equations.  Converting  the  joint  moments  to  cumulants  will  close  this  set  of  equations, 

and  the  numerical  solution  of  such  will  yield  approximations  to  the  elements  of  the  set 

6 

{fCalt...,ae(t)  ■  cii  €  I+,  af  <  3}.  In  order  to  obtain  the  transient  queueing  performance 

*= 1 

measures,  we  note  that 

E[Xi(t)\  =  kei(t),  1  <i<6  (22) 

where  e<  denotes  a  row  vector  whose  ith  element  is  unity  and  all  other  elements  are  zero. 
The  sortie  generation  rate  corresponds  directly  to  the  throughput  of  node  2  of  Figure  2. 
Therefore,  by  applying  Little’s  Law  at  this  node,  the  expected  number  of  sorties  generated 
up  to  time  t,  denoted  by  £[A(t)],  is  given  by 


£[A(t)]  :=  f  fx2ke2(t)dt. 
J  o 
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5  Numerical  Example 


In  this  section,  we  demonstrate  the  implementation  of  the  procedure  described  in 
section  4.  In  particular,  we  set  up  and  numerically  solve  the  set  of  ordinary  differential 
equations  given  by  Eqns.  (20)  and  (21)  for  a  closed  queueing  network  containing  ten  aircraft 
(i.e.,  N  =  10).  As  previously  noted,  the  extension  of  this  model  to  larger  values  of  N 
will  not  increase  the  number  of  differential  equations  in  this  set.  To  further  describe  our 
example  problem,  we  refer  to  the  node  indices  of  Figure  2.  We  assume  that  the  probability 
of  a  pre-sortie  ground  abort  is  pu  —  0.05,  and  we  separately  evaluate  the  model  under  the 
probabilities  p?4  =  0.20, 0.35  of  a  post-sortie  malfunction.  We  assume  that  all  aircraft  are 
initially  grounded  and  waiting  for  taxi  at  node  1;  hence,  kei  (0)  =  10.  All  other  cumulants 
are  set  to  zero  at  time  t  =  0.  All  numerical  results  in  this  section  were  obtained  using  the 
Mathematica®  computing  environment. 

For  demonstrative  purposes,  we  studied  the  expected  workload  and  throughput  of 
the  network  over  the  transient  period  which  is  defined  to  be  the  time  interval  [0,40].  The 
expected  workload  at  each  node  corresponds  directly  to  the  first  order  cumulant  of  X(<),  and 
the  expected  throughput,  i.e.  expected  number  of  sorties  flown  by  time  t,  may  be  calculated 
directly  by  using  Eq.  (23).  Based  on  previous  empirical  testing  in  Matis  [8],  first  order 
cumulant  approximations  under  m  —  3  for  high  traffic  networks  of  this  general  topology  are 
relatively  tight,  i.e.  5%. 

The  output  of  our  numerical  example  for  the  stochastic  routing  probabilities  pu  = 
0.05,  P12  =  0.95,  p24  =  0.20,  and  p2 3  =  0.80  is  displayed  through  Figures  (3)-(9).  In 
particular,  we  have  plotted  the  expected  workload  and  number  of  sorties  flown  as  a  function 
of  time.  It  is  noted  that  each  of  these  measures  approaches  a  steady-state  condition  rapidly 
using  the  particular  parameter  values. 
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Figure  3:  Expected  workload  at  node  1. 


Figure  4:  Expected  workload  at  node  2. 


E  [A  ( t )] 


Figure  9:  Expected  number  of  sorties  by  t. 

The  output  of  our  numerical  example  for  the  stochastic  routing  probabilities  p14  = 
0.05,  P12  =  0.95, ^24  =  0.35,  and  pis  =  0.65  is  displayed  through  Figures  (10)-(17).  It  is 
interesting  to  note  the  dampened  oscillatory  patterns  of  the  measures  E[Xi(t)]  for  i  =  5, 6 
corresponding  to  Figures  (14)  and  (15)  respectively.  The  emergence  of  this  pattern  as  the 
probability  of  a  post-sortie  malfunction  increases  is  likely  due  to  the  network  being  closed 
with  a  limited  number  of  aircraft,  the  large  differences  in  the  parameters  ut(n)  with  n  as 
listed  in  Table  3,  and  the  truncation  of  cumulants  at  m  =  3.  The  overall  expected  holding 
time  of  the  aircraft  E[X\{t)  +  X2(t)]  —  2?[Ai(i)]  +  E[X2(t)\  of  Figure  (16)  does  not  oscillate 
greatly  as  do  the  individual  nodes,  which  is  typically  the  measure  of  practical  interest. 


E[X1(t)]  E[X2(t)] 


Figure  10:  Expected  workload  at  node  1.  Figure  11:  Expected  workload  at  node  2. 
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Figure  12:  Expected  workload  at  node  3. 


Figure  13:  Expected  workload  at  node  4. 


E  [X5(t)l  E[X6(t)] 


6  Conclusions 


This  work  has  considered  the  important  concept  of  operational  readiness  for  an  in¬ 
dividual  air  base,  particularly  when  the  measures  of  operational  readiness  depend  explicitly 
on  time  due  to  theater-level  dynamics.  This  paper  has  extended  a  previous  closed  queueing 
network  model  for  the  aircraft  sortie  generation  process  by  solving  for  vital,  time-variant 
performance  measures  including  the  workload  in  the  system  as  well  as  the  expected  number 
of  sorties  flown  over  a  transient  time  period.  Moreover,  we  generalized  the  maintenance  ser¬ 
vice  time  distribution  by  considering  a  phase-type  (PH)  representation  that  accommodates 
a  transient  analysis  via  the  method  of  cumulant  functions  and  allows  for  inclusion  of  the 
likely  phenomenon  of  blocking  at  the  repair  hangar.  In  particular,  we  assumed  that  a  single 
repair  hangar  accommodates  at  most  two  aircraft.  By  adapting  and  employing  the  cumu¬ 
lant  function  method,  we  obtained  a  computationally  tractable  set  of  ordinary  differential 
equations  that  are  independent  of  the  number  of  aircraft  in  the  system.  Our  approach  to  the 
problem  allows  for  the  computation  of  mean  performance  measures  as  an  explicit  function 
of  time.  Though  the  cumulant-based  procedure  for  transient  analysis  was  tailored  for  this 
queueing  network  model  of  the  aircraft  sortie  generation  process,  it  is  well  suited  for  other 
applications  in  manufacturing,  transportation,  distribution  networks. 
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